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Abstract 

Based  on  the  fully  nonlinear  Boussinesq  equations  in  Cartesian  coordinates,  the 
equations  in  generalized  coordinates  are  derived  to  adapt  computations  to  irregularly- 
shaped  shorelines,  such  as  harbors,  bays  and  tidal  inlets,  and  to  make  computations 
more  efficient  in  large  nearshore  regions.  Contravariant  components  of  velocity  vec¬ 
tors  are  employed  in  the  derivation  instead  of  the  normal  components  in  curvilinear 
coordinates  or  original  components  in  Cartesian  coordinates,  which  greatly  simplifies 
the  equations  in  generalized  curvilinear  coordinates.  A  high-order  finite  difference 
scheme  with  staggered  grids  in  the  image  domain  is  adopted  in  the  numerical  model. 
The  model  is  applied  to  five  examples  involving  irregular  coordinate  systems.  The  re¬ 
sults  of  these  cases  are  in  good  agreement  with  analytical  results,  experimental  data, 
and  the  results  from  the  uniform  grid  model,  which  shows  that  the  model  has  good 
accuracy  and  efficiency  in  dealing  with  the  computations  of  nonlinear  surface  gravity 
waves  in  domains  with  complicated  geometries. 
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1  Introduction 


Boussinesq  models  for  surface  gravity  waves  have  been  proved  to  be  effective  tools 
to  simulate  wave  propagation  in  coastal  regions.  Since  the  introduction  of  the  stan¬ 
dard  Boussinesq  equations  for  variable  water  depth  by  Peregrine  (1967),  numerical 
models  based  on  his  equations  have  been  developed  by  many  researchers  (Goring 
1978;  Abbott  et  al.  1979,  1984;  Elgar  &  Guza  1985;  Liu,  Yoon  &  Kirby  1985;  Rygg 
1988)  and  have  been  shown  to  give  good  predictions  in  comparison  with  field  data 
or  laboratory  data,  when  applied  within  their  range  of  validity.  Recently,  extended 
forms  of  Boussinesq  equations  were  derived  by  Madsen,  Murray  &  Sorensen  (1991) 
and  Nwogu  (1993),  among  others,  to  improve  dispersion  relationships  in  intermediate 
water  depths  and  to  simulate  wave  propagation  from  relatively  deep  to  shallow  wa¬ 
ter.  Wei  &  Kirby  (1995)  developed  a  high-order  numerical  model  based  on  Nwogu’s 
equations  and  provided  additional  validation  tests  of  the  model.  More  recently,  fully 
nonlinear  Boussinesq  equations  were  derived  by  Wei  et  al.  (1995).  The  resulting 
equations  not  only  have  improved  linear  dispersion  properties  in  intermediate  water 
depths  but  they  are  not  limited  to  small  amplitude  waves.  A  time-domain  numerical 
model  based  on  the  equations  was  then  developed  and  verified  against  a  broad  range 
of  experimental  data  (see  Wei  and  Kirby,  1998;  Kennedy  et  al.,  1999  and  Chen  et  al., 
1999a).  Chen  et  al.  (1999b)  applied  the  fully  nonlinear  Boussinesq  model  w-ith  the 
incorporation  of  energy  dissipation  by  wave  breaking  to  investigate  the  fully-coupled 
interaction  of  surface  waves  with  rip  currents  and  the  nearshore  circulation  generated 
by  wave  breaking  on  a  barred  beach  with  a  rip  channel. 

Most  of  the  numerical  Boussinesq  models  are  solved  by  the  finite  difference  method 
using  rectangular  grids.  However,  the  geometric  complexity  of  the  general  coastal 
environment  together  writh  the  rapid  change  in  wavelength  as  ’waves  move  from  deep  to 
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shallow  water  makes  the  use  of  undistorted  grids  in  the  models  somewhat  problematic. 
Typically,  a  model  grid  chosen  to  study  the  entire  range  from  deep  water  to  the 
shoreline  results  in  problems  such  that  waves  are  over-resolved  in  deep  water  and 
under-resolved  in  shallow  water.  To  resolve  a  broad  spectrum  of  wind  waves,  the  time- 
domain  Boussinesq  model  with  a  regular  grid  spacing  may  become  too  expensive  to  be 
used  in  large  nearshore  regions.  Moreover,  complicated  geometries,  such  as  harbors 
and  tidal  inlets,  make  the  computation  very  expensive  because  of  local  resolution 
problems,  and  in  addition,  the  stair-stepping  boundaries  associated  with  rectangular 
grids  may  decrease  the  computational  accuracy.  To  deal  with  such  problems,  irregular 
grid  methods  are  generally  used  through  coordinate  transformations,  such  that  the 
physical  irregularly  shaped  domain  is  transformed  into  a  regular  image  domain  where 
the  finite  difference  computations  are  carried  out. 

There  are  numerous  examples  of  irregular  grid  methods  in  the  study  of  waves. 
Numerical  models  based  on  the  parabolic  approximation  of  the  mild-slope  equation 
for  linear  wave  propagation  in  a  non-orthogonal  coordinate  system  were  developed 
by  Tsay  and  Liu  (1982),  Isobe  (1986),  Liu  and  Boissevain  (1988),  Kaku  and  Kirby 
(1988),  Kirby  (1988),  Kirby  et  al.  (1994).  Dalrymple  et  al.  (1994)  used  spectral 
methods  to  study  forward-propagating  water  waves  in  conformally-mapped  channels. 
An  example  of  the  use  of  a  curvilinear  grid  system  to  solve  Boussinesq-type  equations 
is  given  in  Wang  et  al.  (1992),  who  studied  solitary  wave  scattering  by  a  vertical 
cylinder. 

Irregular  grid  methods  are  also  widely  used  in  numerical  modeling  of  large-scale 
oceanographic  problems  involving  shallow  water  equations.  Non-orthogonal  boundary- 
fitted  grid  models  were  developed  (e.g.,  Sheng  1986;  Chen,  Zheng  and  Zhu,  1999) 
for  modeling  coastal  and  estuarine  processes.  Shi  and  Sun  (1995)  developed  a  self- 
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adaptive  grid  model  for  computing  the  time  dependent  moving  boundary  problem  of 
storm  surge  flooding  using  the  shallow  water  equations. 

Generally,  there  are  three  coordinate  transformation  methods  for  the  transforma¬ 
tion  of  hyperbolic-type  equations  (including  Boussinesq-type  equations  if  only  con¬ 
sidering  leading  order  terms).  The  first  method  is  to  transform  the  independent 
variables  only,  retaining  the  unchanged  primitive  variable  (u,  v,  rj)  (see,  for  example, 
Hauser  et  al. ,  1985,  1986;  Raghunath  et  al.  1987;  Borthwick  and  Barber,  1992).  In 
this  method,  the  transformed  expressions  of  both  equations  and  boundary  conditions 
are  more  complex  than  their  Cartesian  counterparts.  It  is  not  efficient  to  use  the 
method  to  treat  the  fully  nonlinear  Boussinesq  equations  since  the  higher  order  dif¬ 
ferential  terms  in  the  equations  would  make  the  expressions  very  complicated.  A 
typical  example  is  that  uxx  may  be  extended  into  six  terms  in  the  curvilinear  coordi¬ 
nate  and  there  are  even  more  terms  for  higher  order  differentials.  The  second  method 
is  the  adoption  of  the  tangential  velocity  in  curvilinear  coordinates  (e.g.,  Chen,  Zheng 
and  Zhu,  1999)  or  the  covariant  component  of  the  velocity  vector  (Warsi,  1998)  in  the 
coordinate  transformation.  The  resulting  equations  from  these  methods  are  relatively 
simple  compared  with  the  equations  from  the  former  method.  However,  difficulties 
were  found  in  realizing  lateral  slip  boundary  conditions  in  non-orthogonal  curvilin¬ 
ear  coordinates.  The  third  method  is  the  contravariant  component  method  (Warsi, 
1998).  The  contravariant  components  of  the  velocity  vector  can  be  regarded  as  gener¬ 
alized  components  of  velocity  in  the  transformed  image  domain.  Several  advantages  of 
using  the  contravariant  velocity  have  been  recognized  in  the  derivations  of  hyperbolic- 
type  equations,  as  shown  by  Sheng  (1986),  among  others.  Shi  and  Sun  (1995,1997) 
introduced  the  contravariant  components  in  their  coordinate  transformation  of  shal¬ 
low  water  equations  and  easily  obtained  the  kinematical  lateral  boundary  conditions, 
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i.e.,  the  contravariant  components  of  velocity  are  zero  at  boundaries.  Furthermore 
Shi  and  Sun  (1998)  derived  a  new  set  of  equations  for  the  shallow  water  equations  in 
terms  of  contravariant  velocity  and  surface  elevation  in  order  to  solve  the  transformed 
equations  by  using  an  altemating-direction-implicit  scheme. 

In  this  paper,  contravariant  velocity  techniques  are  used  in  the  coordinate  transfor¬ 
mation.  Fully  nonlinear  Boussinesq  equations  in  terms  of  contravariant  components 
of  velocity  vector  at  a  reference  elevation  and  the  surface  elevation  are  derived  in 
generalized  curvilinear  coordinates  based  on  the  fully  nonlinear  Boussinesq  equations 
in  Cartesian  coordinates.  Following  the  work  of  Wei  and  Kirby  (1995)  and  Wei  et 
al.  (1995),  a  fourth-order  Adams-Bashforth- Moulton  predictor-corrector  scheme  is 
employed  in  the  numerical  model  to  perform  the  time  integration.  Unlike  the  spatial 
discretization  in  Wei  et  al.  (1995),  we  use  a  staggered  grid  system  in  the  transformed 
image  domain.  The  first-order  spatial  derivative  terms  are  discretized  to  fourth-order 
accuracy  by  using  standard  five-point  finite-differencing,  and  the  dispersive  terms  are 
discretized  to  second-order  accuracy.  This  ensures  that  the  truncation  error  does  not 
contain  terms  which  are  mathematically  similar  to  the  actual  dispersive  terms. 

The  numerical  model  is  then  applied  to  five  cases  involving  irregular  coordinate 
systems.  The  first  is  wave  evolution  in  a  rectangular  basin  with  a  computational  curvi¬ 
linear  grid.  Consistency  is  found  between  the  curvilinear  grid  model  and  a  uniform 
rectangular  grid  model.  The  second  case  is  wave  shoaling  on  a  sloping  beach.  Here, 
the  element  Courant  number  at  every  grid  point  is  kept  the  same  everywhere  through 
adjustments  of  the  grid  size.  The  example  illustrates  the  gains  in  efficiency  afforded 
by  the  method  in  an  open  coastal  application.  The  third  case  is  wave  propagation  in 
a  circular  channel  in  which  a  curvilinear  grid  is  generated.  The  fourth  case  examines 
diffraction  of  a  solitary  wave  by  a  straight  vertical  wall  at  normal  incidence.  Locally 
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fine  grids  are  generated  in  this  case  around  the  tip  of  the  wall  to  resolve  the  scale  of 
the  wall.  This  case  shows  that  the  present  model  is  capable  of  the  computation  of 
the  nonlinear  wave  propagation  with  good  accuracy.  Finally,  the  model  is  applied  to 
Ponce  de  Leon  Inlet  (Florida,  USA).  Boundary-fitted  grids  with  high  resolution  near 
structures  and  inside  the  inlet  are  generated  and  monochromatic  waves  are  simulated. 

2  Equations  in  Generalized  curvilinear  coordinates 

The  fully  nonlinear  Boussinesq  equations  in  Cartesian  coordinates  (Wei  et  al,  1995) 
are  written  in  terms  of  a  reference  velocity  u  =  (it,  v)  at  a  reference  elevation  za.  The 
mass  conservation  equation  can  be  written  as 

%  +  V-M  =  0,  (1) 

where  M  is  the  depth-integrated  volume  flux  given  by 

M  =  (h  +  V)u+(h  +  V)[£-^(h2-hV  +  n2)}V(V-u) 

+(h  +  rj)[za  +  -  *?)]V[V'(/m)j,  (2) 

in  which  r}  is  the  free  surface  elevation  relative  to  the  still  water  level  and  h  is  the 
still  water  depth.  The  associated  momentum  conservation  equation  is 

nt  4-  (u  •  V)u  4  gVrj  4-  Vx  +  V2  =  0,  (3) 

where  g  is  the  gravitational  acceleration  and  Vx  and  V2  are  the  dispersive  Boussinesq 
terms: 

Vi  =  ^V(V-uO  +  zaV[V-(hnt)}  -  V[^2V-ut  4-  (4) 

V2  =  V{(ztt-?7)(u'  V)[V-(hu)]  +  i(4-7?2)(«-  V)(V-u)} 

+\' V{[V-(hu)  +  t?V-u]2}.  (5) 
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A  coordinate  transformation  is  introduced  in  the  generai  form 


6  =  6(*,v)>  6  =  6 


(6) 


where  (£i,£2)  are  new  independent  variables  in  the  transformed  image  domain.  Re¬ 
ferring  to  Figure  1,  boundaries  Fi,  F2,  F3  and  F4  in  the  physical  domain  (x,  y)  become 
nI;n2,  II3  and  n4  respectively  in  the  image  domain  (£i,£2). 

Instead  of  (u.v)  in  Cartesian  coordinates,  the  contravariant  components  of  the 
velocity  vector  at  za  are  introduced  as  the  dependent  variables  in  the  curvilinear 
coordinate.  In  the  tensor  space,  a  velocity  vector  can  be  expressed  as  two  components 
of  any  covariant  basis,  i.e. 

«  =  tfgh  (?) 


where  u*  are  the  contravariant  components  of  the  covariant  basis  g*.  As  a  simple 
example,  (u,v)  may  be  the  contravariant  components  of  the  covariant  basis  (i,j)  in 
Cartesian  coordinates.  In  the  present  curvilinear  coordinates,  u1  =  (U,  V),  in  which 
{U,  V )  are  the  contravariant  components  of  the  covariant  basis  (gi,g2).  According  to 
the  relationships  among  components  in  different,  basis: 


u 


dxk 


dx 


-u 


(8) 


where  ()?  denotes  the  new  basis.  ( U,V )  can  be  described  by  (u,  v)  in  Cartesian 
coordinates  as: 


U 

V  = 


1 


y/9o 
1 


-  vx6), 


y/9<> 


(~mx +vx&), 


where  g0  is  the  determinant  of  the  metric  tensor  defined  by 


(9) 

(10) 


go 


9n  912 
921  922 


(11) 
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and  where  gij  is  the  covariant  metric: 


dx1'  dx1’ 

9 13  ~  dxi  dx? ' 

(12) 

Using  the  relations: 

dy  ,-56  dx  56 

56  ^9o  dx  ’  <96  V9ody 

dy  _ 56  dx  _ 56 

56  “  ~^9o~dx1  56  “  5y 

(13) 

yields  the  definition  of  (U.  V ): 

dx  dy  dt  ‘ 

(14) 

T/  56  ,  56  dg 2 
^  Udx’rVdy  dt ' 

(15) 

Equations  (14)  and  (15)  indicate  that  (£/,  V)  can  be  regarded  as 

generalized  velocities 

in  the  generalized  coordinates. 

From  (9)  and  (10),  we  can  also  get  the  relationship  between  (U,  V )  and  the  velocity 

components  (un,vn)  normal  to  the  curvilinear  coordinates  (Shi  and  Zheng.  1996): 

Vi?0  Tr 

un  =  v/ - u, 

V§22 

(16) 

_  \fOQ  y 

V9n 

(17) 

Equations  (16)  and  (17)  show  that  (U,  V)  are  stretched  velocity  components  normal 
to  the  curvilinear  coordinates.  The  introduction  of  (U,  V )  makes  it  convenient  to 
obtain  the  lateral  boundary  conditions  as  shown  in  Section  4. 


Equation(l)  and  (2)  are  now  written  in  tensor-invariant  forms: 

ru  +  4=  =  °.  (18> 

^f%OXk 
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(19) 


Mk  =  (*+,)»“-!- (A +  ,)[|-i(A2- A, +  ,2)][-LA(vSKl)!!t 
+(h  +  »?)(*«  +  \ (h  -  n)}[-j=^(Vfohu‘)]it. 

The  tensor-invariant  forms  of  the  momentum  equations  (3),  (4)  and  (5)  are: 


(20) 


(21) 


(22) 

where  k,l,m  —  1  and  2,  (u1,  u2)  —  (U.V),  (rrhir2)  =  (Ci- £2);  is  the  partial 
derivative,  ()*  represents  the  covariant  spatial  derivative  (See  Appendix)  while  lk 
represents  the  contravariant  spatial  derivative  defined  as 

f,k  =  | J9 *•  (23) 

3  Numerical  method 

Numerical  analyses  (Wei  and  Kirby  1995)  of  finite  differencing  for  Boussinesq  equa¬ 
tions  show  that  it  is  necessary  to  adopt  high-order  schemes  in  either  space  or  time  in 
Boussinesq  equations  because  the  truncation  errors  of  a  second-order  approximation 
may  contaminate  the  real  dispersive  terms  in  the  equations.  It  also  should  be  noted 


duk 

dt 


f  gr}lk  +  ttVj  +  Vf  -r  V*  =  0, 


tln,k  »  T/& 


v? 


'  1,2  a  (VToul)  + 


l2y^  dxl 


y/9odzl 


V-2 


1 ,  2 
+2(Za 


rf)ul 


d 


dxl  \/(fo  dx7 


+kl^=-i;(^hu‘)  +  JL±(y/g;u‘)n\\ 


2  y/godx1 


■s/godx1 
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that  the  finite  differencing  on  irregular  grids  can  lead  to  a  loss  of  accuracy  if  dis¬ 
cretizations  of  the  untransformed  equations  are  performed  on  irregular  grids.  In  the 
present  paper,  we  can  avoid  the  problem  by  discretizing  the  transformed  equations 
in  the  image  domain  with  regular  grids. 

By  using  the  tensor  formula  presented  in  the  Appendix,  the  system  of  equations  is 
rewritten  in  a  form  that  makes  the  application  of  the  difference  procedure  convenient 
as  that  of  Wei  et  al.  (1995).  The  mass  conservation  equations  (18)  and  (19)  may  be 
expressed  as: 

r)t  =  E{rl,U,V),  (24) 

where 

E  =  +  »)E%,  +  k Mb  +  n)\ %} 

y/90 

+  >?)  +  -  V 2)}^z(DU)Sl 

*  V9o 

— 7={(ci h2(h  +  i f)  +  - r}(h 2  - 

y/9o  b  y/9o 

+{a2h(h  +  n)  -  rf(h  +  rj)\-~z{DHU)^}^ 

+-]={{alh*{h  +  17)  +  \r}{h*  -  n ?)}M;(DU)s2 
VPo  6  v9o 

+\a2h(h  +  ,)  -  \n (ft  +  n)\^L{DHV)b}(l 

+[02h(h  +  ,)  -  \n(h  +  ,)J £jL(DHU)b }«„  (25) 

in  which, 

DU  =  -±=(VToUh  +  (26) 

V3o  v9o 

DHU  =  -^=(V^hD%  +  4=(VSf hV )6.  (27) 

V3o  V5o 
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The  momentum  equation  in  direction  can  be  written  as 


where 


in  which 


Ut  =  F(rl,U,V)^\F1{y%  +  {F2(V,V)}t 


+Ft(n,  ut,  vt)  +  Ut,  Vt)  +  ft(u,  u,  V)  +  Fein,  V,  V),  (28) 

V  =  U  +  A26,%-i=(v^kk  +  hb2f-[^=(VmhU)(,h,,  (29) 

9o  V5o  5o  v^o 

F  =  — j-(p22%  -  5,12^2)  “  +  ^^2 

go 

+DlnU 2  +  2D\2UV  +  D\ 2V2),  (30) 

5c  V.9o 

-A62^[-|=(VSAV)fe]5..  (31) 

5o  V5o 

fi(tl,V)  =  fc2f>,f^(v»'k  +  -4=(^)&]Si 

5o  Vw  V5o 

+M2^[-^=(v^M06  +  -UveJMOfck,  (32) 

50  V^O  yj  9 0 

Fs(v,Ut,Vt)  =  ^[^(DUTJ+^DHOT)]*,  (33) 

0o  2 

Ft(n,ut,vt)  =  -^[^(PUD+d(I>BUT)]6,  (34) 

fife  £7,  K)  =  -^[(*,  -  n)U(DHV%  +  (za  -  n)V(DHU)b 
go 

+i(4  -  >!2MSt%  +  i(4  -  V2)V(DU)b k 

-i^[((PFf7)+IJ(DU))2]£„  (35) 

^  <?o 

Fs(v,U,V)  =  ^-[(za  -  rf)U{DHU)^  +  (za  -  V)V(DHU)(2 
9o 

+1(4  -  n*)U(BU)t,  + |(4  -  4)v(oc/-)&]& 

+i^(((i?OT)  +  -)(I>tO)2)ei,  (36) 

Ot/r  =  -2=(V9^<)«.  +  -^=(VSSVi)&.  (37) 

V5o  V5o 
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DSUT  =  -L(^hUt)Si  +  -^={^rnhVt)b.  (38) 

v50  y/9o 

In  equation  (28),  we  introduce  F,  F1}  ...  .  F6  to  represent  separately  the  terms  with 
different  properties.  F  includes  the  pressure  gradient  terms  and  convective  terms; 
F\  and  F2  are  the  linear  dispersive  terms,  in  which  F>  is  the  additional  term  due  to 
the  non-orthogonality  of  the  coordinate;  F3  and  F4  are  the  nonlinear  dispersive  terms 
with  time  derivatives  while  F4  is  the  additional  term  due  to  the  non-orthogonality; 
F=  and  F6  are  the  nonlinear  dispersive  terms  with  spatial  derivatives  while  F6  is  the 
term  due  to  the  non-orthogonality.  Djk  is  the  Christoffel  symbol.  The  constants 
Gi;  Go,  61, 62  are  defined  as 


ax  =  01  /2  -  1/6,  a2  =  /?  +  1/2,  b1  =  f,b2  =  8  (39) 

where  8  =  zajh,  0  =  —0.531  in  the  present  paper. 

Similarly,  the  momentum  equation  in  £2  direction  can  be  written  as 

Vt  =  G{r},U,V)  +  [Gl{U)}tF{G2{U,V)]t 

FGs(V:  Ut,  V))  +  G4(v,  Ut ,  Vt)  +  Gs(v,  u,  V )  +  G6(v,  U,  V),  ’  (40) 


where, 

v  =  V  +  h%f-{^(^)(,\h+h^{-±(s/FehV)b}(1,  (41) 

9o  y/9o  9o  v5o 

G  =  --j-(-gi2ri^  +  giiVti)  ~ 

9o 

+D2nU 2  +  2D%lUV  +  F|2V'2),  (42) 

Gi(to  « 

9o  \/9o 

-W^/UvSSM'ls,  ]&,  (43) 

9o  \/9o 

<h(V,V)  =  +  /=(vW0&k, 

Po  vao  VFo 
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+M2^[^=(v^W)£l  +  -^:(V5bAV) &j£„  (44) 

2o  v^o  v^o 

G3(<7,Di,VJ)  =  &-[£{DUT)+n(DmjTj\b,  (45) 

^0  £ 

G4(v,Ut,Vt )  =  [^(I>£7T)  +  rj(DHUT)]€l,  (46) 

a5(r?;  (/,  V)  =  ~^[(aQ  -  n)U(DHV%  +  (*a  -  J?)y(Dro)$, 

<?o 

+\(4  -  V2)U(DL%  +  i(4  -  ^(DR)^, 

-—[((OJT^+'lPO-))2)^,  (47) 

G6(I),  R  V)  =  ^[(z„  -  i j)U(DHU)b  +  (za  -  v)V(DHL% 

ffo 

+1(4  -  V2)U(Dl%  +  i(4  -  ti‘)V(DU)b]b 
-i^[(  (DHU)  +  n(DU))%.  (48) 


G,  Gi,  ...  ,  Gg  have  the  similar  meaning  as  F ,  F1;  ...  ,  F6  as  described  above. 

The  arrangement  of  cross-differentiated  and  nonlinear  time-derivative  terms  on 
the  right  hand  side  of  equations  (28)  and  (40)  makes  the  resulting  set  of  left-hand 
sides  purely  tridiagonal. 

A  staggered  grid  in  —  £2  plane  is  employed  as  shown  in  Figure  2,  where  the 
crosses  denote  rj  -  points  at  which  rj  is  computed,  the  circles  denote  (/-points  at 
which  U  and  U  are  computed  and  the  squares  denote  V-points  at  which  V  and  V 
are  computed.  The  separate  labeling  for  rj.  U  and  V  in  this  scheme  is  convenient  for 
the  implementation  of  boundary  conditions.  The  new  co-ordinates  (£x,£2)  are  taken 
as  integer  grid  positions,  £1  =  1, 2. ....  m,  £2  =  1,2,  Following  the  work  of  Wei 
and  Kirby  (1995),  we  discretize  the  first-order  spatial  derivative  terms  to  fourth-order 
accuracy  by  using  standard  five-point  finite-differencing,  leading  to  truncation  errors 
of  0(A£4).  The  dispersive  terms  themselves  are  finite-differenced  only  to  second-order 
accuracy,  leading  to  error  terms  of  0(A£2)  relative  to  the  actual  dispersive  terms. 


13 


The  fourth-order  Adams-Bashforth-Moulton  predictor-corrector  scheme  is  em¬ 
ployed  to  perform  time  updating.  A  sequence  of  time  instants  are  defined  by  t  =  pAt. 
Level  p  refers  to  information  at  the  present,  known  time  level.  The  predictor  step  is 
the  third-order  explicit  Adams-Bashforth  scheme,  given  by 


;23(£)?,  -  16(B)?;1  +  5(B)?;2], 


W 


%  +  ^[23 (F)fc  -  16(F)?;1  +  5(F)?;2), 
%  +  |£[23(G%  -  16(G')?;'  +  5(G')?;2], 


where 

F'  =  F  -fi  (Fi)t  +  (ih)*  +  F$  +  F4  +  F*,  +  F6,  (52) 

G'  —  G  -f  4-  (<?2 )t  +  O3  +  G4.  +  G$  +  G$.  (53) 

In  equations  (50)  and  (51),  {Fi)t,  {F2)t,  (Gi)t  and  (G2)t  involves  time  derivatives  and 
the  calculation  of  F3,  F4,  G?J  and  <?4  at  a  certain  time  level  requires  the  corresponding- 
values  of  Ut  and  Vt.  They  can  be  evaluated  by 


= 

= 

MS2  = 


2^[S <  ~  4<- 1  +  <J2]  +  0(At% 
±-Hj-w??]  +  0(At% 

"2 k[Zwfj2  ~ 4wfjl + <i] + 0(At% 


(54) 

(55) 

(56) 


where  w  represents  U,  V,  F\:F2,  G\  or  G2. 

When  Ufj1  and  V?j~l  are  obtained  from  equations  (50)  and  (51).  the  contravariant 
velocity  (U,  V)  at  the  new  time  level  can  be  solved  by  a  system  of  tridiagonal  matrix 
equation,  using  elimination  method. 

After  the  predicted  values  of  (rj,  U,  V)?]1  are  evaluated,  the  corresponding  quan¬ 
tities  of  (E,  F',  G%j  at  time  levels  (p  + 1),  (p),  (p  - 1),  (p  -  2)  are  obtained.  Then  we 
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use  the  fourth-order  Adams-Moulfon  corrector  method: 


vr» 


At 


=  ^  PWff  +  19(E)!,  -  5(B)?-1  +  (£)?;2], 


IJV. 

uh3 


24 

At 

24 

At 


[9 (-ns1  +  19(*%  -  5(F)fJx  +  (F)fJ 


r^P-2] 


%  +  24  PCTfT  +  19(G')?4  -  S(G')?;1  +  (G')f-2). 


(57) 

(58) 

(59) 


Similar  to  the  treatment  in  the  predictor  stage.  Ut,  Vt,  (Fi)t,  (F2)t,  ( Gi)t  and  {G-2)t  are 
evaluated  in  the  following  manner: 


II 

r~< 

eL  (UaF  18<i  +  9<7  ”  2“0  +  0(A‘3)’ 

(60) 

Mfc  = 

t(H'  +  H  ~  KJ1  +  +  o(Ai>), 

(61) 

M  CJ1  = 

~e + 3v*F 1  -  K* + <f  I  +  0(Ai!), 

(62) 

11 

c* 

1  **3 

-g|j(ll«<72  -  18<J*  +  -  2<f)  +  0(A«3). 

(63) 

The  corrector  step  is  iterated  until  the  error  between  two  successive  results  reaches 
a  required\limit.  The  error  is  computed  for  each  of  the  three  dependent  variables 
rj,  U,  V  and  is  defined  as  in  Wei  et  al.  (1995): 


A  /  = 


g ij  i/f. 


■JH-1 


l/f, 


(64) 


where  /  denotes  77.  £7  or  V'  and  ()*  denotes  the  previous  results.  The  corrector  step 
is  iterated  if  any  of  the  A /  exceeds  10-4. 

Initial  testing  of  the  present  staggered-grid  scheme  has  shown  it  to  be  a  consider¬ 
able  improvement  over  the  scheme  of  Wei  and  Kirby  (1995),  in  terms  of  short  wave 
length  noise  generation.  A  more  extensive  exploration  of  the  behavior  of  the  two 
schemes  is  underway  and  will  be  reported  separately. 
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4  Boundary  conditions 


For  computations  of  wave  propagation  in  domains  with  complicated  boundaries,  ap¬ 
propriate  boundary  conditions  have  to  be  specified  in  the  numerical  model.  Generally, 
for  a  perfectly  reflected  vertical  wall,  the  horizontal  volume  flux  normal  to  the  wall  is 
zero.  Wei  and  Kirby  (1995)  derived  a  set  of  boundary  conditions  satisfied  by  normal 
velocity  at  zQ,  tangential  velocity  and  surface  elevation  by  considering  only  the  lead¬ 
ing  order  terms  in  the  mass  conservation  equation,  i.e.,  for  the  case  of  a  vertical  wall 
parallel  to  x  axis  ,  the  boundary  conditions  are 

v  =  0;  rjy  =  0;  uy  =  0.  (65) 

They  then  applied  five-point  off-center  finite  difference  to  the  equations  above.  These 
boundary  conditions  may  be  adopted  in  curvilinear  coordinates  if  we  use  the  con- 
travariant  component  (U,  V)  instead  of  {u,v).  That  is 


Cl 

II 

0 

at  IL,n2, 

(66) 

11 

0 

at  n3,n4, 

(67) 

77!1  =  0;  K!1 

=  0  at  nx,n2, 

(68) 

Tjx?  =  0;  Ul2 

=  0  at  XI3 ,  II4 . 

(69) 

The  boundary  conditions  (66)  and  (67)  can  be  easily  obtained  in  the  staggered  grids 
because  U  points  are  located  on  Ey,  II2  and  V  points  on  II3,  II4.  The  symmetric  or 
anti-symmetric  conditions  are  used  to  apply  the  boundary  conditions,  i.e.,  U  values 
are  anti-symmetric  to  Hi  and  n2,  V  values  are  anti-symmetric  to  II3  and  II4  according 
to  (66)  and  (67),  while  U  and  r)  values  are  symmetric  to  II3  and  IL-,  V  and  r\  values 
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are  symmetric  to  III  and  II2  according  to  (68)  and  (69).  The  symmetric  or  anti¬ 
symmetric  conditions  also  make  it  easy  to  discretize  equations  on  boundaries  or  in 
the  vicinity  of  boundaries,  especially  for  high  order  finite  differences  with  five-point 
schemes. 

'The  sponge  layer  boundary  condition  and  the  wave  generating  boundary  condition 
are  also  implemented  in  this  paper  following  Wei  and  Kirby  (1995). 

5  Examples 

5.1  Wave  evolution  in  a  rectangular  basin 

As  a  simple  yet  efficient  test  case  (Wei  and  Kirby,  1995),  the  evolution  of  waves 
in  a  rectangular  basin  was  calculated  by  using  both  the  uniform  rectangular  grid 
model  and  the  curvilinear  grid  model.  Though  there  are  no  corresponding  nonlinear 
analytical  solutions  or  experimental  data  to  compare  with,  the  comparison  between 
the  results  from  the  two  models  can  show  a  consistency  of  the  solutions  on  curvilinear 
grids  and  rectangular  grids. 

The  basin  dimensions  are  20m  x  20m,  and  the  water  depth  is  0.5m  constant  over 
the  basin.  The  initial  condition  is  provided  by  a  motionless  Gaussian  hump  of  water 
with  its  center  located  at  the  center  of  the  basin  (xc,yc): 

7j(x, y,t  =  0)  =  Ho exp{-7[(a:  -  xcf  +  {y-  yc)2]},  (70) 

u(x,  y,  t  =  0)  =  0,  (71) 

v(x,  y,  t  =  0)  =  0,  (72) 

where  Ho  is  the  initial  height  of  the  hump,  7  is  the  shape  coefficient,  and  (xc,  yc) 
is  the  coordinate  at  the  center  of  the  domain.  We  chose  Hq  —  0.2 m,  7  =  0.4,  and 
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Zc  =  Vc  =  10  m. 

In  the  rectangular  grid  model,  grid  sizes  are  set  to  be  identical  everywhere  as  Ax  = 
Ay  —  0.1  m.  For  the  curvilinear  grid  model,  a  grid  generation  method  (Brackbill, 
1982)  can  be  used  to  generate  the  curvilinear  grid.  As  a  test,  we  generate  a  curvilinear 
grid  by  weighting  the  center  point  as  shown  in  Figure  3  in  which  the  maximum  grid 
size  is  0.15m  near  boundaries,  while  the  minimum  grid  size  is  0.045m  at  the  center 
of  the  domain.  The  total  number  of  grid  points  is  200  x  200,  as  the  same  as  that  of 
the  regular  grid  model. 

Figure  4  shows  water  surface  contour  in  the  physical  domain  and  the  contour 
obtained  from  the  rectangular  model.  There  is  very-  small  difference  between  the 
dashed  line  and  the  solid  line  in  Figure  4.  The  small  differences  are  caused  by  the 
interpolation  errors  and  by  the  different  resolutions  of  the  initial  Gaussian  hump. 
Numerical  experiments  show  that  the  difference  decreases  with  a  reduction  of  grid 
sizes  in  both  models. 

5.2  Monochromatic  waves  shoaling  on  a  slope 

As  a  second  example  using  the  coordinate-stretching  approach,  we  study  periodic 
wave  propagation  over  a  sloping  beach  where  there  is  a  contraction  of  the  wavelength 
and  a  resulting  increase  in  required  resolution  as  waves  shoal  towards  shore.  The 
previous  uniform  grid  model  generally  makes  the  computation  expensive  because  fine 
grids  with  constant  grid  spacing  are  always  adopted  in  order  to  resolve  short  waves 
in  shallow  water,  and  further,  the  time  step  needs  to  be  small  in  order  to  have  an 
appropriate  Courant  number  in  deep  water.  The  problem  can  be  solved  by  using 
an  irregular  grid  in  which  the  grid  spacing  is  adjusted  according  to  the  water  depth 
so  that  the  Courant  number  in  each  element  is  constant  in  the  whole  domain.  The 
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idea  has  been  used  by  Kashiyama  and  Okada  (1992)  in  shallow  water  flow  analysis. 
To  show  the  efficiency  of  the  present  model  in  dealing  with  wave  propagation  over 
a  sloping  beach,  we  present  a  1-D  model  with  the  stretched  grid  size  in  the  present 
case. 

Figure  5  shows  a  mild-slope  beach,  where  a  constant  depth  of  8m  on  the  left 
connects  to  a  constant  slope  on  the  right.  Waves  with  a  period  of  4s  are  generated 
by  a  wavemaker  located  on  the  left  side.  Two  sponge  layers  are  placed  at  both  ends 
of  the  domain  to  absorb  wave  energy.  The  uniform  grid  model  and  the  stretched  grid 
model  are  then  used  respectively  in  this  case. 

In  the  uniform  grid  model,  a  constant  grid  size  of  0.4m  is  adopted,  giving  a  total 
of  2500  grid  points  in  the  computation  domain.  According  to  linear  theory,  the 
corresponding  wavelength  in  the  deepest  water  is  45.2m,  while  shortest  wavelength 
in  the  shallow  water  is  18.4m.  Thus  there  are  46  grid  points  per  wavelength  in  the 
shallowest  region  on  the  right  side  and  113  grid  points  per  wavelength  in  the  deepest 
water  region  on  the  left  side.  The  choice  of  an  adequate  resolution  in  shallow  water 
thus  leads  to  an  over-resolution  of  the  wave  form  in  deep  water.  In  addition,  for 
a  given  time  step,  the  Courant  number  increases  with  the  increase  of  water  depth. 
Therefore,  a  time  step  has  to  be  selected  according  to  the  Courant  number  in  deep 
water.  Here,  the  time  step  in  the  uniform  grid  model  is  chosen  to  be  0.025s  which 
leads  to  the  Courant  number  of  0.19  in  the  shallowest  region  and  0.471  in  the  deepest 
region. 

In  the  stretched  grid  model,  gradually  varying  grid  sizes  are  chosen  from  1m  in  the 
deepest  region  to  0.40m  in  the  shallowest  region  according  to  the  following  coordinate 
transformation: 

f  (73) 


19 


where  C  is  the  wave  celerity  from  the  linearized  Boussinesq  equations  and  A  is  a 
coefficient,  A  =  7.52  in  this  case.  The  Courant  number  is  constant  in  the  whole 
computational  domain.  The  total  grid  number  decreases  to  1420,  43%  less  than  that 
of  the  uniform  grid  model.  Even  so,  the  resolution  of  the  model  is  not  reduced  since 
there  are  45  grid  points  per  wavelength  both  in  deep  water  and  in  shallow'  water, 
comparing  to  the  46  points  per  wavelength  in  shallow'  water  in  the  uniform  grid 
model.  The  time  step  for  the  stretched  grid  model  can  be  much  larger  than  that  for 
uniform  grid  model  because  of  the  larger  grid  spacing  in  deeper  water.  For  instance, 
the  time  step  in  this  case  is  chosen  to  be  0.06s,  which  is  2.4  times  larger  than  that 
in  uniform  grid  model  and  gives  a  Courant  number  of  0.452.  The  decrease  in  grid 
number  and  increase  in  time  step  make  the  stretched  grid  model  much  more  efficient 
in  comparison  with  the  uniform  grid  model.  In  such  a  case,  more  than  a  four  time 
speed-up  may  usually  be  expected  in  view'  of  the  grid  numbers  and  time  steps  used  in 
the  two  different  models.  The  actual  computational  time  in  the  stretched  grid  model 
is  3.1  times  faster  than  the  uniform  grid  model,  which  is  a  little  slower  than  expected 
because  the  stretched  grid  model  needs  a  few  more  iterations  in  the  calculation.  In 
the  computation  of  wave  propagation  on  a  natural  beach  with  a  relatively  mild  slope 
offshore  and  a  steep  slope  close  to  the  shore  line,  the  efficiency  will  be  more  obvious. 

The  surface  elevations  obtained  from  the  stretched  grid  model  are  presented  in 
the  image  domain  as  shown  in  Figure  6.  It  is  clear  that  the  wavelength  in  the  image 
plane  is  almost  constant,  which  illustrates  that  the  model  provides  the  same  resolution 
when  waves  propagate  from  deep  water  to  shallow  water.  The  surface  elevations  in  the 
physical  domain  are  shown  in  figure  7  in  comparison  with  the  results  from  the  uniform 
grid  model.  Figures  8  and  9  respectively  display  the  comparisons  of  wave  height  and 
wave  number  betwreen  the  two  models.  The  wave  numbers  here  are  obtained  from 
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the  wave  elevation  distributions  calculated  in  the  two  models.  The  comparisons  show 
that  the  stretched  grid  model  results  are  in  good  agreement  with  results  from  the 
uniform  grid  model. 

5.3  Waves  in  a  circular  channel 

The  present  example  is  a  constant  depth  channel  with  vertical  sidewalls  laid  out  in 
a  circular  planform.  Dalrymple  et  al  (1994)  used  spectral  methods  with  coordinate- 
transformed  equations  to  analytically  study  linear  wave  propagation.  In  their  study, 
three  cases  with  different  widths  of  channel,  namely,  narrow,  wider  and  very  wide 
channels,  w-ere  carried  out.  The  case  of  the  very  wide  channel  presents  a  more  com¬ 
plicated  pattern  of  waves  with  diffraction  and  strong  reflection.  Thus,  this  case  is 
chosen  for  study  in  this  paper,  and  comparison  is  made  between  the  numerical  results 
and  analytic  solution. 

Let  ri  and  7*2  be  the  inner  and  outer  radius  of  the  channel  respectively,  with 
7*1  =  75m,  r2  =  200m  in  this  case.  The  depth  of  the  channel  is  4m.  The  coordinate 
transformation  can  be  described  as: 


£2  =  han  *(f)  (75) 

7T  X 


The  grid  mesh  in  the  physical  domain  is  shown  in  Figure  10.  The  grid  spacing  in  the 
radial  direction  is  constant  at  1m,  while  a  constant  angular  grid  is  used  along  the 
channel  length,  resulting  in  a  maximum  tangential  grid  size  of  1.26m  near  the  outer 
wall  and  a  minimum  of  0.47m  near  the  inner  wall.  The  waves  propagate  primarily 
counter-clockwise  from  the  mouth  of  the  channel.  A  wavemaker  is  located  internally  at 
the  eastern  end  of  the  channel  and  two  sponge  layers  are  placed  in  straight  channels 
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extending  from  the  eastern  and  western  end  of  the  circular  channel  (not  shown  in 
Figure  10).  The  wavemaker  produces  linear  waves  with  a  period  of  4.s  and  a  very 
small  amplitude  for  comparison  with  the  linear  analytic  solution.  The  time  step  is 
chosen  as  0.05s  in  the  present  case. 

To  illustrate  the  process  of  wave  propagation  in  the  channel,  Figures  11-13  show 
the  transient  propagation  of  a  wave  train  into  the  channel.  It  is  shown  that  the  waves 
initially  propagate  in  a  straight  line,  but  as  the  channel  bends,  the  waves  start  to 
diffract  around  the  bend  and  simultaneously  run  into  the  curving  channel  sidewall 
and  are  reflected  around  the  bend.  The  present  method  allows  for  transient,  wave 
propagation  while  the  spectral  method  given  by  Dalrymple  et  al.  can  only  describe  a 
steady-state  linear  solution.  Figures  14  and  15  depicted  the  comparisons  of  the  water 
surface  variation  along  outer  wall  and  inner  wall  between  the  analytic  solution  and 
numerical  solution  after  a  periodic  steady  state  has  been  achieved.  Good  agreements 
are  found  in  the  comparisons. 

5.4  Diffraction  of  a  solitary  wave  by  a  straight  vertical  wall 

An  experiment  on  the  diffraction  of  a  solitary  wave  by  a  straight  thin  wall  was  car¬ 
ried  out  by  Perroud  (1957)  in  a  wave  tank.  A  diffracting  wall  of  0.8cm  thickness  is 
placed  in  the  tank  to  diffract  the  solitary  wave  at  normal  incidence.  Figure  16  shows 
the  experiment  layout,  where  the  measurements  were  made  at  a  point  (p,  6)  in  the 
defined  polar  coordinate.  The  tests  were  run  at  a  constant  water  depth  of  6.1cm  and 
various  ratios,  ranging  from  0.27  to  0.58,  of  the  incident  wave  height  to  water  depth. 
Measurements  of  wave  heights  and  patterns  of  the  diffracted  waves  were  performed 
in  different  directions  8  from  the  diffracting  wall  and  the  central  axis  at  the  distance 
P- 
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To  simulate  the  diffraction  of  waves  in  the  experiment,  it  is  necessary  to  use 
stretched  grid  model  to  resolve  the  thin  diffraction  wall.  Thus  a  stretched  rectangular 
grid  mesh  is  generated  in  the  computation  domain  with  linear  variation  of  grid  sizes 
along  the  axis  direction.  The  finest  grids  are  obtained  around  the  wall  with  the 
smallest  grid  size  of  0.8cm,  which  is  exactly  the  thickness  of  the  wall.  The  largest 
grid  size  in  the  whole  domain  is  2.54cm  (1m)  at  the  two  sides  of  the  domain.  The 
grid  size  in  the  normal  direction  of  the  axis  is  identical  as  2.54cm.  The  location  of 
initial  solitary  wave  is  set  at  76.2cm  eastward  from  the  wall. 

The  diffraction  and  reflection  process  after  an  incident  solitary  wave  impacts  on 
the  wall  is  shown  in  a  time  sequence  of  contour  plots  of  the  surface  elevation  for  the 
case  of  ao/h  =  0.42  (Figure  17).  When  the  -wave  impinges  on  the  wall  (t  =  0.8s), 
the  wave  on  the  upward  side  runs  up  and  reflects  back.  The  initial  development  of 
the  diffracted  wave  can  be  clearly  seen  near  the  tip  of  the  wall.  Following  the  wave 
impact,  the  early  stages  of  propagation  of  the  diffracted  solitary  wave  and  the  back- 
scattered  reflected  wave  are  then  shown  at  t  =  0.9 s.  At  t  =  1.4s,  the  initial  primary 
wave  separates  into  a  forward  transmitted  diffracted  wave  and  a  backward  scattered 
reflected  wave.  The  newly  evolved  secondary  baekscattered  and  forward-scattered 
•waves  generated  from  the  tip  of  the  -wall  propagated  outward  and  follow  the  leading 
reflected  and  diffracted  waves,  respectively. 

The  numerical  results  and  measurement  data  obtained  from  four  solitary  waves 
with  different  wave  heights  are  compared  in  Figure  18  (a)- (h).  The  diffraction  co¬ 
efficient  k,  defined  as  the  ratio  of  the  diffracted  wave  height  to  the  incident  wave 
height  (a/ao),  is  plotted  against  the  distance  p  for  different  values  of  the  angle  9. 
The  close  agreement  proves  that  the  present  model  is  capable  of  simulating  nonlinear 
wave  propagation. 
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5.5  Application  to  Ponce  de  Leon  Inlet 

To  demonstrate  the  practicality  of  the  present  model,  we  use  the  mode!  to  simulate  the 
propagation  of  monochromatic  waves  in  Ponce  de  Leon  Inlet,  Florida,  USA.  Smith 
and  Harkins  (1997)  used  numerical  models  to  estimate  wave  transformation  over 
Ponce  de  Leon  Inlet  and  comparisons  were  made  between  numerical  results  and  mea¬ 
surement  data  from  the  physical  model  The  geometry  is  shown  in  Figure  19,  which 
consists  of  a  coastline,  a  jetty,  an  inlet  leading  to  the  Halifax  and  Indian  Rivers,  and 
a  complex  bathymetry.  A  boundary-fitted  grid  is  generated  as  shown  in  Figure  20.  In 
order  to  resolve  structures  and  short  waves  in  shallow  water,  finer  resolution  is  used 
near  the  jetty,  coastal  lines  and  the  inlet.  At  the  offshore  boundary  of  the  domain, 
monochromatic  waves  are  generated  -with  a  period  of  15s  and  a  small  amplitude.  To 
avoid  wave  reflections  by  the  coastal  boundaries  and  wave  breaking  in  shallow  water, 
sponge  layers  are  placed  in  shallow  water  area  along  the  coastlines.  Figure  21  shows 
a  snapshot  of  surface  elevations,  showing  wave  reflection  on  the  upward  side  of  the 
jetty,  wave  diffraction  on  the  leeward  side,  refractive  wave  focusing  in  the  area  to  the 
right  of  the  inlet  mouth,  and  standing  waves  inside  the  inlet.  Though  the  numerical 
results  have  not  been  compared  with  measurement  data  as  there  is  no  consideration  of 
energy  dissipation  by  wave  breaking  and  wave  runup  on  sloping  beaches  yet,  this  case 
study  illustrates  that  the  present  model  has  a  potential  prospect  for  computations  in 
complicated  domains,  such  as  in  harbors  and  tidal  inlets. 

6  Conclusions 

Based  on  the  fully  nonlinear  Boussinesq  equations  derived  by  Wei  et  at  (1995),  the 
equations  in  generalized  curvilinear  coordinates  are  derived  by  using  contravariant 
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velocity  method.  Then  the  numerical  model  is  developed  with  a  high-order  finite 
difference  scheme  in  staggered  grids.  To  test  the  model,  five  examples  involving 
curvilinear  or  stretched  coordinate  systems  are  applied.  The  computation  of  wave 
evolution  in  a  rectangular  basin  with  a  curvilinear  grid  indicates  that  the  model  is 
consistent  with  the  uniform  grid  model.  In  the  case  of  wave  propagation  on  a  sloping 
beach,  the  same  element  Courant  number  is  obtained  in  every  grid  points  by  adjusting 
the  grid  size,  so  that  the  resolutions  of  waves  are  the  same  both  in  shallow  water  and 
deep  water.  The  computational  efficiency  is  shown  to  be  greatly  improved  by  the 
new  model.  Wave  propagation  in  a  circular  channel  is  simulated  by  employing  the 
present  model  with  the  boundary-fitted  grid.  Good  agreement  is  found  between  the 
numerical  results  and  analytic  solution.  Then  the  model  is  used  to  simulate  the 
diffraction  of  a  solitary  wave  by  a  straight  vertical  wall  at  normal  incidence.  The 
comparison  between  numerical  results  and  measurements  shows  that  the  model  has 
good  accuracy  in  dealing  with  the  computation  of  nonlinear  wave  propagation  with 
complex  lateral  boundaries.  Finally,  the  model  is  applied  to  Ponce  de  Leon  Inlet.  In 
this  case,  monochromatic  waves  are  simulated  in  the  complex-shaped  domain  with  a 
real  bathymetry. 

For  practical  application  in  complicated  domains,  the  present  model  needs  to  be 
further  improved  with  the  incorporation  of  energy  dissipation  by  w-ave  breaking  and 
w-ave  runup  on  sloping  beaches  and  structures.  The  development  will  be  reported  on 
in  the  near  future,  in  conjunction  with  a  more  complete  investigation  of  Ponce  de 
Leon  Inlet. 
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Appendix:  Spatial  derivative  in  tensor  forms 

The  tensor-invariant  forms  of  controlling  equations  (18)  -  (22)  can  be  expanded 
into  equations  (24)  -  (48  )  by  using  the  following  tensor  formulas. 

The  gradient  of  a  scalar  /  can  be  written  as: 

v/  =  fg‘  P6) 

where  gi  is  the  contravariant  basis.  According  to  the  relationship  between  two  differ¬ 
ent  basis: 

g‘  «  9iSS j  (77) 

in  which  is  the  contravariant  metric,  the  gradient  can  be  expressed  on  the  covariant 
basis  gi  as 

vf  =  =  f-Z‘  (78) 

In  the  present  paper,  for  example,  since  the  contravariant  velocity  is  introduced  as 
the  dependent  variable,  the  equations  have  to  be  expended  on  the  covariant  basis  g,. 
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The  pressure  gradient  term  may  be  expressed  as 
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where  £0  is  the  determinant  of  the  metric  tensor.  In  equation(79),  we  employed  the 
following  formula: 

W‘  =  S{  (80) 


where  8?  is  the  Kronecker  delta. 

The  divergence  of  a  vector  u  can  be  written  as 

1  dy/lfru1 
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In  the  present  paper,  V-u  can  be  expanded  as 

1  d 


(81) 


V-u  = 
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where  U  and  V  are  the  contravariant  components  of  velocity  vector. 

To  a  vector  u  in  a  two-dimensional  space,  the  covariant  spatial  derivative  Uj  can 
be  expressed  in  terms  of  the  contravariant  component  u1  of  u  as  following 
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where  Dljk  is  the  ChristoSel  symbol.  Then  the  convective  terms  in  the  paper  can  be 
expanded  as  follows: 
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Figure  5:  Schematic  of  the  wave  flume. 
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Figure  6:  The  stretched  grid  model:  surface  elevations  in  the  image  domain 


Figure  7:  Surface  elevation  comparison  between  the  uniform  grid  model  (solid  line) 


and  the  stretched  grid  model  (dashed  line)  in  the  physical  domain. 
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Figure  8:  Comparison  of  wave  heights  between  uniform  grid  model  (solid  line)  and 
stretched  grid  model  (dashed  line). 


Figure  10:  Calculation  grid  in  the  circular  channel. 
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Figure  11:  Wave  propagation  in  the  circular  channel (t=  50  seconds). 
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Figure  12:  Wave  propagation  in  the  circular  channel(t=  100  seconds) 
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Figure  13:  Wave  propagation  in  the  circular  channel{t=  200  seconds) 
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Figure  14:  Comparison  of  the  surface  variation  along  outer  wall  between  the  analytic 
solution  (solid  line)  and  numerical  solution  (dashed  line). 


Figure  15:  Comparison  of  the  surface  variation  along  inner  wall  between  the  analytic 
solution  (solid  line)  and  numerical  solution(dashed  line). 
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Figure  20:  Computational  grid  in  Ponce  de  Leon  Inlet 
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Figure  21:  A  snapshot  of  wave  surface  elevations  in  Ponce  de  Leon  !nlet(T  =  15s) 
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